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I present time-dependent analytical solutions of quasi-steady shocks with cooling, where quasi-steady shocks are objects composed of 
truncated steady-state models of shocks at any intermediate time. I compare these solutions to simulations with a hydrodynamical 
code and finally discuss quasi-steady shocks as approximations to time-dependent shocks. Large departure of both the adiabatic 
1 and steady-state approximations from the quasi-steady solution emphasise the importance of the cooling history in determining the 
trajectory of a shock. 
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1 Introduction 



Analytic solutions of well defined problems are often used as ben chmark tests for hydrodynamical codes. 
However, t he most widely used tests such as the Sedov blast wave (|Sedovlll993l ;i or the Sod shock tube test 



are almost all adiabatic or with Mach numbers of order unity. By contrast, gas in astrophysical 



contexts can be subject to strong cooling and Mach numbers of order a few hundreds are common place 
in protostellar jets. There is hence a lack of analytical solutions with cooling or at high Mach numbers. I 
present here time-dependent analytical solutions of quasi-steady shocks with cooling for arbitrarily high 
Mach numbers. 

i iLesaffre et all (|2004al ) and ILesaffre et al. I (j2004bl ) studied the temporal evolution of molecular shocks 
^""^ ' and found that they are most of the time in a quasi-steady state, ie: an intermediate time snapshot is 
composed of truncated steady state models. They showed that if the solution to the steady state problem 
Jlf^ , is known for a range of parameters, it is possible to compute the time-dependent evolution of quasi-steady 

^ ■ shocks by just solving an implicit ordinary differential equation (ODE). Since the steady state problem is 

O itself an ODE, it is very easy to numerically compute the temporal evolution of quasi-steady shocks. 
^ However, a given cooling function of the gas does not necessarily lead to an analytical steady state 
solution. Furthermore, even when it does, the implicit ODE which drives the shock trajectory does not 
necessarily have an analytical solution itself. In this paper, we tackle the problem in its middle part : 
d ' we assume a functional form for the steady state solutions (section 2). We then show how to recover the 
underlying cooling function that yields these steady states (section 3) . Finally, we exhibit a case where the 
shock trajectory has an analytical solution (section 4) and we compare it to numerical simulations (section 
5). Results are discussed in section 6 and conclusions drawn in section 7. 



2 Method 



Consider the following experimental set up: we throw gas with a given pressure, density and supersonic 
speed fo on a wall. We assume a perfect equation of state with adiabatic index 7. We assume the net 
rate of cooling of the gas is a function A{p,p) where p and p are the local density and pressure of the 
gas. The gas is shocked upon hitting the wall, heated up by viscous friction and an adiabatic shock front 
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develops that soon detaches from the wall. Behind this front, the gas progressively cools down towards a 
new thermal equilibrium state and a relaxation layer unrolls. 

All physical quantities are normalised using the entrance values of pressure and density, so that the 
sound speed of the unshocked gas is cq = ^/7. The time and length scales will be later specified by the 
cooling time scale (see section I3.ip . 

Consider now the set of all possible stationary states for the velocity profile in the frame of the shock 
front. A given entrance speed uq in the shock front provides the velocity u at a given distance y behind 
the shock front: 

u = f{y,uo). (1) 

The adiabatic jump conditions for an infinitely thin (or steady) shock enforce /(0,uo) = Ua{uo) where 

7-1 27 1 

■Ua(wo) = — — r^^o + — - — • (2) 
7 + 1 7 + 1 Mo 

Unfortunately, that / is a simple algebraic function of y and uq does not necessarily imply an algebraic 
form for A{p,p). It is in fact more appropriate to express y in terms of u and uq. We hence rather write 
([1]) in the following manner: 

y = g{u,uo) (3) 

with the condition 

g{ua,uo) = 0. (4) 
S ection 13.11 detai l s how t o recover A{p,p) from g{u,uo 



Lesaffre et al. ( 2004bl ) provide the ODE for the evolution of the distance from the shock to the wall 



with respect to time t if the shock is quasi-steady at all times: 

r = f{r,r + vo) (5) 

where a dot denotes a derivative with respect to time. This equation can also be expressed with the 
function g: 

r = g{f,r + vo). (6) 

In sectional show how one can integrate equation Q and give an analytical expression for a simple form 
of g. In sectional compare this solution to a time-dependent numerical simulation. 

3 Cooling function 
3.1 General procedure 

Let us write the equation of steady-state hydrodynamics in the frame of the shock with entrance parameters 
{p,p,u) = (1, l,no): 

pu = uo, (7) 
p + pu"^ = l + ul (8) 
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and 

|W:^P + ^P-^)]=A (9) 

with the boundary condition u = Ua{uo) at y = 0. 

One can solve the equations ([7]) and ([8]) for p and p and use the relations into ^ which becomes: 

d 1 

A{u,y,uo) = ^[^^(— ^-(1 - uuo + ul) + -uuq)] (10) 
ay 7 — 1 2 

Expansion of the derivative with respect to y provides: 

A{u,u',uo)= ^ (I + ul)u' - ^ ~^ ] uquu' (11) 
7 — 1 7 ~ 1 

where u' = du/dy. 

By taking the derivative of equation ([3]) we easily extract u' in terms of u and uq: 

u' = l/^{u,uo). (12) 
combined with p!T]) provides A(u, uq). ([7]) and ([8]) finahy give A{p,p). 



3.2 First application 

I now illustrate this method with a simple function g{u, uq). In a typical radiative shock in the interstellar 
medium, the post-shock gas finally cools down to nearly the same temperature as the pre-shock gas. To 
mimic this effect, we need a cooling function such that its thermal equilibrium is the isothermal state. In 
other words, A{p,p) = implies p = p. This is equivalent to asking that the final steady velocity of any 
shock verifies the isothermal jump condition u = Ui where Ui{uo) = 1/uq: 

lim 7r^(u, no) = — oo. (13) 
To verify both conditions (jU and (|13p we take the simple form: 

g{u, no) = P -— (14) 

u-Ui (no ) 

where /3 > determines the strength of the cooling. Setting a length scale allows to assume [5=1. The 
above procedure (section [3.ip yields: 

A(p.,) = -lf (l+7) |l+p(7P-7-l)](p-rt l\ 

i^V',-if(p-iY'(vp--i)~j7W^)] 

This cooling function is displayed on figured) In addition to the temperature T = p/ p = 1 solution, the 
thermal equilibrium state A{p,p) = is also realised when the factor [1 + p (7 p — 7 — 1)] is set to zero. 
However, this state is practically never achieved in the relaxation layer of a shock as it happens for densities 
p< 1 + 1/7 and p is always greater than the adiabatic compression factor Cq = (7 + l)/(7 — 1 + 27nQ ^) ~ 
1 + 2/(7 — 1) for strong shocks. In the high temperature limit A scales like which is reminiscent of 
the collisional coupling between gas and dust for low dust temperatures. But in the high density limit, 
A ~ p^2T^ which yields a rather unphysical scaling on the density. In the next subsection, I show how 
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Figure 1. Contour plot of the cooling function A defined by equation I I15I I for 13=1 and 7 = 5/3 with respect to temperature T = p/p 

and density p. Levels for the contours are (solid line) and respectively -1, -10, -10^ and, -10^ for dotted, dashed, dash-dotted and 
dash-dot-dot-dotted lines. Overlayed (long-dash) are three typical shock thermal profiles at age t = 100 (from bottom to top vq = 10, 

21 and 46). 

to improve the physical relevance of the cooling function, at the loss of the analytical integrability of the 
trajectory. 



3.3 Second application 

In this section, I briefly illustrate how one can obtain semi-analytic approximation of shocks for any kind 
of cooling function. I start with a given cooling function Aq and compute an analytical approximation to 
the steady state function g{u,UQ). I then recover the underlying cooling function Ai for this approximate 
steady state and check how Ai and Aq differ. 
A very common form for the cooling due to a collisionally excited (unsaturated) line is 



13 



(16) 



where Tq is the temperature of the transition and (3 scales the strength of the cooling. I use (3 = 1 for 
simplification (it amounts to specify a length scale without loss of generality). 
If we apply the procedure of section 13.11 backwards we have to integrate 



dy 
du 



u 



[7(1 + ul) - (1 + 7)nno] 



(7 - 1) < 



expl 



TqUq 



M (1 + tig — UUq) 



(17) 



to find the equation for the stationary velocity profile. This equation does not have an analytical solution. 
But there are many approximations to the right hand side that will allow to treat the problem. 

For example, we can simplify ()17p by using the strong shock approximation uq » 1 along with the high 
compression approximation u » uq: 



dy 
du 



7M 

7^ 



exp[ 



UUq 



(18) 
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Figure 2. Contour plot of the cooling functions Aq and Ai for f3 = 1, Tq = 1000 and 7 = 5/3 with respect to temperature T = p/p and 
density p. Levels for the contours are (solid line, only for Ai) and -lO'^, -1, -10~^, -10~ and -10~® from top to bottom for dotted 

(Ao) and dashed (Ai) lines. 



This equation still does not have an analytical solution but we can add a term to the factor of the 
exponential to get the derivative of the simple function G{u,uq) = ii^exp[To/(niio)]. Hence we finally take 

— = (1-4-^) exp[ ] (19) 

dn J 7 — 1 ^0 

which will be a good approximation provided that To >> Ug. This yields the simple form G{u,uo) — 
G{ua,UQ) for the function g{u,uo) with which equation Q can then be integrated numerically and tabu- 
lated to get a fast access to the shock trajectory. 

To check that the above approximations did not alter too much the underlying cooling function, we can 
apply the procedure 13.11 to the simplified equation ([19]). This provides: 

7(p- 1)[4 - 4p + To(p - 1)] p-1 

Figure E] compares contour plots of both cooling functions Ag and Ai (for Tq = 1000 and 7 = 5/3). 
It can be seen that despite the crude approximations we made, Aq and Ai are still close to one another 
for a very large range of parameters (for 1 < T << Tq and p >> 1 both expressions asymptotically 
converge). However, thermal equilibrium solutions (solid lines in figure [2]) appear for Ai when none existed 
for Aq. Also, the range of applicable entrance velocities is restricted to conditions such that the maximum 
temperature in the shock is low compared to Tq (because we made the Tq » Uq approximation). 

This is nevertheless a good illustration of this method which can in theory be applied to any cooling 
function. Indeed, one can in principle use the fact that u/uq is bounded by a constant number strictly 
lower than 1 to uniformly approach equation (llTh with polynomials (or any other functional form simple 
to integrate). One then recovers an analytical expression for g(u,UQ) arbitrary close to the exact solution. 
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4 Shock trajectory 
4.1 Exact solution 

Implicit ODEs like ([6]) are in principle straightforward to integrate numerically. It is however much harder 
to find an analytically integrable form for these equations. Such a solution nevertheless exists for the simple 
but physically relevant example (I14p . 
Let us use into ([6]) with /? = 1 to obtain 

^_ 27 + (7-l)z;g + (7-3)rz;o ^^i) 
(7 + 1) (r^ ^rvQ- 1) 

The solution of (j2ip for r yields only one positive root r = h{r): 



h{r) = (l-^)^o-{l + ^>^r + s{r) 
^ ' 2(7 + l)r ^ ' 



where 



s(r) = \/a + br + cr'^ (23) 

with 

a=(7-3)2^;2, (24) 

6 = 2(7 + 1)[47 + (7 + 1)t;2] (25) 

and 

c = (^ + 1)2(4 + ^2)_ (26) 

We are hence able to express the age t of the shock as a function of its position r which provides the 
trajectory: 

i(0 = r 7^ (27) 



Hx) 



We now write l/h{x) as a sum of integrable terms: 



' -^0+ / + .+^-^^+ + (28) 



h{x) 2 4(7 + l)s(a;) 2s{x) 2{z + x) 2{z + x) s{x 
with 

d=:^vo{3 + vl) , (29) 
7 + 1 

e = 47 + 5(7 - l)vl + (7 - 1) ^0 (30) 
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and 

27 , 7-1 2 



+ ^< (31) 



7+1 7+1 

We then integrate expression psp to obtain the trajectory in its final analytical form 
, , vo s(r) — ^/a e , / b-\-2^/ac \ 



2 2(7 + 1) 2^ ^\b + 2cr + 2^/^s{r)) 2 ^V^ + ^ 

(7 + l)d^ / (z + r)(2a-bz + 2^^(-z)) \ 
2s(-z) ^ V^[(2a-6z) + (6-2cz)r + 2s(r)s(-z)]y 

If r[t) is wanted, t(r) can be numerically inverted by a Newton method of zero finding. This is done 
easily since the derivative t'{r) = l/h{r) is provided analytically. 

4.2 High Mach number approximation 

We can recover a more simple but approximate trajectory if we make two assumptions : 

• in the limit of high Mach numbers {vq/cq — > oo), the adiabatic compression factor becomes a constant: 

Ua{uo) ~ '^^—^uq. (33) 
7 + 1 

• at late times, the compression factor is nearly the isothermal compression factor and r ~ I/uq. Hence 
for high Mach numbers = t^o + — ^^0 and we use: 

ni(no) ~ l/vQ (34) 

With both these approximations, (fT^ with j3 = 1 becomes 

I N (l-7)^^o + (7 + l)^i 
g[u,uo) = vq— — — — — , 35 

(7 + 1) {uvo - 1) 



the shock front velocity is: 



Kr) = (7-l)^g + (7+l)^ (3,) 



and the resulting trajectory is 



^7-1 2 ^ ^, J (7-1)^^, 



t{r) =vor + vo i^v'^ - — -) log \' .\ \ . (37) 

7 + 1 7 + 1 V(7 - l)'^o + (7 + 1)?^/ 

For early times (small r) the strong adiabatic shock trajectory is recovered: 

t{r) ^ (38) 
7 - 1 t;o 

For very late times (very large r) the isothermal shock trajectory is reached asymptotically: 

t{r) ~ vq r. (39) 
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5 Numerical simulation 
5.1 Numerical method 



I compute here the time evolution of a radiative shock with coolin g function (|15p thank s to a ID hy- 



drodynamical code. This code makes use of a moving grid algorithm ( Dorfi &: Drury 19871 ) which greatly 
helps to resol ve the adiabatic froii t while keeping the total number of zones fixed to 100. The mesh driving 
function (see Dorfi &: Drurvl 198?! ) is designed to resolve the temperature gradients. 



The advection scheme is upwind, Donnor-Cell. The time integration is implicit fully non-linear with an 
implicitation parameter of 0.55 as a compromise between stability and accuracy. The time-step control 
keeps the sum of the absolute values of the variations of all variables lower than 0.5. In practice, the 
maximum variation of individual variables at each time-step is lower than 1%. 

I use a viscous pressure of the form: 

= -pc,v'(Ax/10)2 + Z2max(-^, 0) (40) 
3 ox 



where Cg = y/TP/p is the local sound speed. Ax is the local grid spacing and I = 10~^ is a prescribed 
dissipation length. 

To avoid numerical difficulties due to the form of A{p,p) when p or p are close to 1, we set A to zero 
when p — 1 OT p — 1 are lower than 10~^. The time normalisation is such that (3=1. 

The entrance parameters are set to {p,p,vo) = (1,1,10) with 7 = 5/3 and the evolution is computed 
until a stationary state is nearly reached. 



5.2 Trajectory 

The position of the shock at each time-step is computed as the position of the maximum of the ratio Pv/p 
along the simulation box. I compare this trajectory to the analytical expression (j32|) on figure [3l At a given 
position, the relative difference on the ages of the shock is maximum at the very beginning, when the shock 
front is being formed and the position is still no more than a few dissipative lengths. For times greater 
than 5 x 10~^, the relative error is less than 8%, with a secondary maximum at r ~ 1. An estimate for this 
error is given in section 16.11 Note that both the isothermal (j39p and the adiabatic ([38]) approximations are 
wrong by an order of magnitude at this point. The adiabatic approximation is accurate to a 20% level only 
for ages lower than 0.1. Afterward the effects of cooling slow down the shock and the adiabatic solution 
overestimate the position by large. The isothermal approximation is valid up to a 20% level only for times 
greater than 3000. This is because it does not take into account the period of early times when the shock 
is moving swiftly and since the isothermal shock moves at a slow pace it takes time to recover the delay. In 
other words, the cooling history of the shock does make a difference to its position. Since the ratio between 
the adiabatic speed and the isothermal speed scales like Vq this situation will be even worse for stronger 
shocks. By contrast, the error estimate given in section [6T] suggests that the quasi-steady approximation 
is equally good for stronger shocks. 

The approximate trajectory for high Mach number quasi-steady shocks ()37p is also shown. It is already 
a very good approximation even for the relatively low Mach number (vq/cq = 7.75) I use. Indeed in section 
14.21 1 neglected only terms of order greater than 2 in the inverse of the Mach number. 



5.3 Snapshots 

I output the results of the simulations at a few selected time-steps. For each of these snapshots, I determine 
the position r of the shock with Pv/p as in the previous subsection 15. 2[ I then compute the velocity of the 
quasi-steady shock front r = h{r) thanks to (j22p . uq = vq + r gives the entrance velocity in the frame of 
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Figure 3. Trajectory of the shock age vs position in the simulation (solid) compared to the analytic expression for a quasi-steady shock 
(dashed) and the high Mach number approximation (dotted). Also shown (solid curve in upper pannel) is the relative error between the 
solid and dashed curves as well as the estimate ^ /r for this error given in section 16.11 (dotted curve in upper pannel) . 



the steady shock. I now recover the relation u = f{y,UQ) thanks to equations ^ and (|14p : 

_ -27 + (l-7)^g + (7 + l)j/ 

(7 + 1)^0(^-1) • ^ ' 

The temperature (T = p/p) profile can finally be retrieved from this velocity profile thanks to the relations 
dZD and dS]). 

I compare the quasi-steady state solution to the results of the numerical simulation on figure |4l The 
gas is flowing from the right onto the wall on the left. At early times, when the shock front is still close 
to the wall, the temperature at the wall is higher in the simulation than in the quasi-steady shock. This 
is mainly due to the wall heating effect, which decreases at later times when the cooling function has a 
stronger influence on the temperature. The decrease of the maximum temperature in the shock is due to 
the decrease of the relative entrance velocity of the gas in the adiabatic shock front (see figure H]) . Note 
the high resolution provided by the moving mesh at the adiabatic shock front. For later times at the end 
of the relaxation layer, the temperature decreases toward its final value of T = 1 at equilibrium. 

The dotted curves in figure H] are the high Mach number solutions (see section r4.2p . They already are 
very close to the exact solutions as the approximation is of order 2. 



6 Discussion 



6.1 Time- dependent shocks 

The differences between the quasi-steady state and the numerical solution described in section [5] both 
come from the numerical errors in the scheme and from the fact that quasi-steady shocks are only an 
approximation to time-dependent shocks. I give here an estimate on the difference between quasi-steady 
shocks and shocks. 

I now write the inviscid equations of time-dependent hydrodynamics that a time-dependent shock would 
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Figure 4. Temperature (T = p/p) profiles in the hydrodynamical simulation (diamonds) compared to the analytical solution for a 
quasi-steady shock (solid red line) and its high Mach number approximation (dotted green line). 



set to zero: 

dp d{pv 



^ dv dv 1 dp 
^ dt dx p dx 



and 

^ 1 dp I dipv) dv ^ , . ,,,, 

where p{x,t), v{x,t) and p{x,t) are the density, velocity and pressure fields in the frame of the wall. 
If we now use the quasi-steady solutions to express these equations, we find 



dp. 



= TT^ ^' > (45) 
duo 



E, = {^-l)r (46) 
duo 



and 

E = ^ (47) 
' (7 - 1) duo ^ ^ 

where Ps{y,uo), Us{y,uo) and Ps{y,uo) are the solutions of the steady state equations in the frame of the 
shock. 



2008 3:4 Geophysical and Astrophysical Fluid Dynamics lcsafTre4 



Time- dependent analytic solutions of quasi-steady shocks with cooling 11 

Hence quasi-steady shocks are in general only approximations to time-dependent shocks. The quasi- 
steady state approximation amounts to neglecting the acceleration r of the shock front. Note that a 
maximum departure of the trajectory from the numerical simulation occurs around r ~ 1 when the shock 
switches from adiabatic velocities to isothermal velocities (see figure [3]), ie: when accelerations are likely to 
be the highest. A rough estimate for the relative error on the position is hence rt^/r which overestimates 
the error by more than about a factor 3 (see figure [3]) . Interestingly, this estimate does not depend on the 
shock velocity for strong shocks. 

It is also interesting to note from equations (I45p - (l47p that a high dependence of the steady state on the 
entrance velocit y un will cause departu res of time-dependent shocks from the quasi-steady state. In this 



context, refer to lLesaffre et al. I (|2004bl ) who found that the quasi-steady state was violated for marginally 
dissociative shocks. For this type of shocks, the entrance velocity is indeed close to the critical velocity at 
which the major cooling agent is dissociated and a small variation of the entrance velocity can strongly 
affect the post-shock. 



6.2 Entrance parameters 

In section 3, I used a normalisation based on the entrance values for the density and pressure: this must 
not hide the fact that the cooling function (jlSp implicitly depends on these parameters. Hence for a given 
cooling function, I computed analytical solutions only for a fixed set of entrance density and pressure. At 
this point, there is no reason why other sets of parameters should provide integrable solutions. 



7 Summary and future prospects 

I described a general way of obtaining time-dependent analytical solutions to quasi-steady shocks with 
cooling. I applied this method to a physically sensible example and compared the resulting quasi-steady 
shock to a time-dependent hydrodynamic simulation. I also provided a more simple high Mach number 
approximation to the exact solution. I showed that even though quasi-steady shocks are not strictly 
speaking time-dependent shocks, they are a good approximation to time-dependent shocks. In particular, 
more simple approximations such as the adiabatic approximation or the steady-state approximation badly 
fail to reproduce the behaviour of the shock for a large range of times. This is because the cooling history 
of the shock is essential in determining the position of the shock. 

I wish to emphasise also that the method described in section 13.11 allows to recover the underlying 
cooling function for any set of steady state solutions. The associated quasi-steady shocks can then be 
easily computed by solving the ODE ([6]). As demonstrated in subsection 13.31 one can fit any given cooling 
function with analytical forms for the functions g{u, uq). This hence provides a potentially powerful method 
to quickly compute the evolution of any quasi-steady shock with cooling. 

Analytical solutions of quasi-steady shocks can be used as a basis to study properties of time-dependent 
shocks. Linear analysis around the quasi-steady solution may provide insight for the time-dependent be- 
haviour of shocks. They can also help to address under what conditions shocks tend toward the quasi-steady 
state. 

Furthermore, this work represents a first step towards exact solutions or better approximations to time- 
dependent problems of non-adiabatic hydrodynamics. In the future this might lead to new algorithms for 
dissipative hydrodynamical simulations. Finally, note t hat similar procedure s can be applied to shocks 
with chemistry and magnetic fields (using the results of Lesaffre et al. 2004bl ) . 
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